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Abstract. Populations of competing biological species exhibit a fascinating interplay 
between the nonlinear dynamics of evolutionary selection forces and random 
fluctuations arising from the stochastic nature of the interactions. The processes 
leading to extinction of species, whose understanding is a key component in the study 
of evolution and biodiversity, are influenced by both of these factors. 

In this paper, we investigate a class of stochastic population dynamics models 
based on generalized Lotka-Volterra systems. In the case of neutral stability of the 
underlying deterministic model, the impact of intrinsic noise on the survival of species 
is dramatic: it destroys coexistence of interacting species on a time scale proportional 
to the population size. We introduce a new method based on stochastic averaging 
which allows one to understand this extinction process quantitatively by reduction to 
a lower-dimensional effective dynamics. This is performed analytically for two highly 
symmetrical models and can be generalized numerically to more complex situations. 
The extinction probability distributions and other quantities of interest we obtain show 
excellent agreement with simulations. 
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1. Introduction 

Interactions between biological species are known to lead to very diverse and intricate 
behaviour of a population. This includes, just to give a few examples, coexistence of a 
surprisingly high number of competing species in the same ecological niche [1] , oscillating 
population cycles [2] and chaos [3] . The question which (if any) of the species in a web 
of interactions survive, and for how long, is thus very nontrivial but central for the 
understanding of evolution and biodiversity [4]. 

A classical and long-established model for the interaction of species is the Lotka- 
Volterra model [5, 6]. Since its introduction, it has also been successfully applied in 
many different contexts outside of population dynamics: among others, neural networks 
[7] , game theory [8] and physiology [9] . This model attempts to describe the interaction 
between S species through a set of coupled ordinary differential equations of the form 



The abundance of each species is given by a continuous, real-valued variable Xi with 
i — 1, S . hi are constant source terms describing the growth (or decline) of each 
species in the absence of the others, and Aij is a constant matrix modelling the 
interactions between the species. Within this model, survival or extinction of species 
is purely deterministic: Any fixed initial condition determines unambiguously which, if 
any, of the species survive. Technically, the main underpinning for this is the stability or 
instability of certain stationary solutions of the differential equations (1). Some rather 
precise criteria for determining the persistence of species directly from the vector 6, the 
matrix A and the initial conditions have been obtained in literature [10, 11]. 

However, in a real biological situation, the population is made up from a large but 
still finite number of individuals. Hence, the abundances of each species can only change 
in discrete steps and not continuously. Furthermore, the interactions between them, as 
well as birth and death processes, have - to some extent - a stochastic nature. All these 
features cannot be modelled by the deterministic equations (1). In fact, such effects 
of finite system size and fiuctuations due to some intrinsic randomness (or, likewise, 
due to external noise) have recently been recognized to be very important for extinction 
processes, especially in the case when the deterministic solutions exhibit neutral stability 
(see, for example. [12, 13, 14, 15, 16, 17, 18, 19, 20, 21] and many others). 

In this paper, we propose a new method based on the idea of stochastic averaging 
X which allows one to gain a quantitative understanding of the stochastic extinction 
process in the case when the deterministic limit of the model is neutrally stable. 

In section 2 we will formulate a stochastic model of population dynamics based on 
a graph of interactions between species, whose dynamics in the absence of noise reduces 
to a Lotka- Volterra model of the form (1). 

I The idea of stochastic averaging was first introduced by Khasminskii in [22]. Later on, it was 
rigorously justified in [23] for two-dimensional systems possessing a conservation law. So far, however, 
it has not gained a lot of popularity in physical literature. 




(1) 
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Figure 1. Subgraph of a complicated interaction graph 



In sections 3 and 4 we will treat two pedagogical examples of such models, the 
three-species and four-spccics systems with cyclic dominance. Their deterministic 
dynamics will be shown to be neutrally stable, i.e. to lead to perpetual coexistence 
of all species with periodically oscillating abundances. However, we will see that taking 
into account fluctuations due to the stochastic nature of the interactions introduces 
a finite mean extinction time proportional to the population size. Using stochastic 
averaging, we will characterize the extinction process by an effective stochastic process 
in the deterministically conserved quantities. This will allow us to obtain quantitative 
results on extinction times and their dependence on the initial conditions. 

The generalization to more complex models will be discussed in section 5. 

2. Stochastic Lotka- Volterra Models 

Let us consider a population with 5* species, and interactions between them defined by 
a graph G — (V, E) with vertices V — {Xi, Xs} and (directed) edges 



An edge from Xi to Xj is denoted by an arrow in figure 1, and is taken to mean that the 
species Xi dominates over Xj. We allow at most one edge between each pair of species. 
From this graph, we can now write down a set of reaction equations implementing the 
interactions of the species. For every edge {Xi,Xj) e E we formulate an interaction 
between X^ and Xj in the formalism of chemical reaction equations: 



where the reaction rate kij is the probability for this reaction to occur per (infinitesimal) 
time unit dt and per possible pair of individuals. 

Note that the model described by reactions (2) provides an individual-based, 
discrete description of the interactions, and includes stochastic fiuctuations since the 
reaction rates are interpreted probabihtistically. 



E = {ix,,x,)\ij^r, i,j = i,...,s}. 
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Since the reactions (2) keep the total number of individuals fixed, we can assume 
a constant system size N . The system state is then described by the S'-tuple of species 
counts n = (ni, n^), with rii G 0, ...,A^ and X^^i = We define the 5" "basis 
vectors" ei,...,es by Cj = (0, 0, 1, 0, 0) where the 1 is on the i-th position. The 
reactions (2) can then be translated into a master equation giving the evolution of the 
occupation probabilities Pfi{t) for each state n: 

dtPn{t) = Yl ^ij [i'^i - + ^)Pn-ei+ej{t) - niUjPnit)] . (3) 

{Xi,Xj)&E 

For biological applications, one is mostly interested in large populations, i.e. in 
the limit of large N. The relative abundances of each species, = can then be 
assumed to be real- valued variables in the interval [0; 1]. Using a standard Kramers- 
Moyal expansion [24], the master equation (3) can then be approximated by a Fokker- 
Planck equation in the intensive variables x^: 

d,P{{x,},t) = - Y:^^[a^P{{x,},t)] + ^ d,dAB^jP{{xk},t)]. (4) 

i=l j,J=l 

The conservation of the total population size N gives rise to the normalization condition 
— 1. The drift and noise terms in (4) are given by: 

s 

ctj Xj ^ ] Aj^jXjj (5) 

^ ^ i J2k=i\^ik\xiXk for i = j 
'•^ \ -\Aij\xiXj for iy^j 

The entries of the interaction matrix A are: 

h, if {x,,Xj)eE 

-hj if {Xj,Xi)eE . (7) 
otherwise 

As is well known [25], the Fokker-Planck equation (4) can be reformulated as a set 
of stochastic differential equations §, or Langevin equations: 

1 

dxi — aidt-\ — -= ^ dj dWj . (8) 

C is a matrix satisfying CC^ = B, with B defined by (6). Certainly, this does not fix C 
uniquely, but the precise choice has no influence on the stochastic process [25]. Wj are 
independent Wiener processes, or Brownian motions, with zero mean and unit variance. 

From equation (8) we see that the deterministic, noiseless hmit of the general model 
(2) is given by the following set of coupled ordinary differential equations (so-called rate 
equations) : 



Aij 




dtXi = Xi ^ij^jj ■ (9) 

Throughout this paper, we take ah stochastic differential equations to be in the Ito interpretation 
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This can now be immediately recognized as a generalized Lotka- Volterra model of the 
form (1). The source terms vanish, i.e. hi — Q ||, and the interaction matrix A is given 

by (7). 

The noise terms in (8), proportional to encapsulate the fluctuations due to the 
discreteness of the individuals and the stochastic nature of the reactions (2). 

In total, we have derived a procedure allowing us to obtain a stochastic model (in 
terms of Fokker-Planck or Langevin equations) of a system with a large population of 
individuals from a general interaction graph of the species. This prescription is certainly 
not unique, however, it has the nice property that the deterministic dynamics of the 
resulting model is in one-to-one correspondence with a Lotka- Volterra model whose 
interaction matrix is the adjacency matrix of the original graph. 

Considering the widespread use and the importance of Lotka- Volterra models, it 
seems worthwhile to study models of the form (8) and the effects of stochasticity in 
them. 

In general, the deterministic rate equations (9) possess extinction fixed points 
(where some species are extinct, i.e. there are some j with Xj = 0) and coexistence fixed 
points (where all species are present, > for all i). Prom (9) we see immediately that 
the coexistence fixed points form a linear subspace of the phase space of the system which 
is given by the kernel of the matrix A. When the full stochastic model (8) is considered, 
fluctuations cause the system to touch an extinction hyperplane where Xi = for some 
i sooner or later. Since a species which has died out cannot be re-introduced (this is 
apparent from the reaction equations (2)), this means that in the stochastic system, 
extinction always occurs eventually. 

However, the time scale on which this process occurs can vary greatly. A 
classification of the possible scenarios, characterized by the scaling of the mean 
extinction time T with the population size N was developed in [26] and [15] : 

(i) Stable coexistence for T oc e^, occurring when the deterministic dynamics has a 
stable attractor in the coexistence region. Here, extinction is driven by rare large 
deviations and hence the extinction times for large populations are extremely long. 

(ii) Unstable coexistence for T oc log N. This occurs when the fiow of the deterministic 
dynamics approaches one of the extinction hyperplanes for large times, and weak 
fiuctuations are already sufficient to make one of the species go extinct. 

(iii) Neutrally stable coexistence for a power- law dependence, T oc A^'''. This occurs when 
the deterministic dynamics possesses a family of neutrally stable, closed orbits, 
corresponding to the existence of a conservation law. 

For simple models, these criteria correspond to (linearly) stable, unstable, or neutrally 
stable coexistence fixed points in the rate equations (2). 

II Intuitively, this happens since all reactions in (2) have exactly two reactants. A model with non- 
vanishing source terms can be built by considering reactions of the form Xj — >■ and Xi 2Xi in 
addition to (2) 
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Figure 2. Interaction graph for the rock-paper-scissors game 



Observe that the effects of stochasticity are most dramatic in a neutrally stable 
model: while the deterministic dynamics predicts perpetual coexistence far away from 
the extinction planes, inclusion of fluctuations introduces a finite mean extinction time 
which only scales linearly with the population size. In the following two sections, we 
will now analyze the stochastic extinction process for two pedagogical example models 
of this kind. 

3. Cyclic Three-Species Model: The Rock-Paper-Scissors Game 

3.1. Introduction 

The first example that shall be considered in detail is a three-species model with cyclic 
dominance, whose interaction graph is shown in figure 2. One of the most popular 
areas where such cyclic, intransitive relationships between three entities arise is in game 
theory as a so-called rock-paper-scissors game [8]. In a more biological context, they 
have been observed between strains of E. coli bacteria [27] and between lizard morphs 
[28]. Another rather different application is to forest fire models [29], where the three 
states trees, fire and ash obey a similar relatioship. 

The reaction equations for this model, according to the general treatment in 
section 2, read: 



A + B ^ A + A. 



B+C^B+B 



C + A^C + C. 



(10) 



The interaction matrix is, accordingly: 

/ 1-1 

{A,)= -1 1 
V 1 -1 . 
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Figure 3. Phase space of the three-species rock-paper-scissors game. Grey curves: 
closed deterministic orbits, given by p = const.. Blue (or dark grey) curve: example 
of a stochastic trajectory, obtained for a system size of = 200. 



In order to simplify the calculations, we have set all reaction rates to be equal in 
(11). By rescaling time we can then set them to 1 without loss of generality. According 
to the general treatment in section 2, the stochastic model in the large- limit is then 
described by the Fokker-Planck equation (4) (or, equivalently, the Langevin equations 
(8)) with S = 3. The drift term (5) and the noise term (6) evaluate explicitely to the 
following expressions: 



Note that here and in the following we shall use the variable names a, b, c and xi, 
X2, X3 interchangeably. A qualitative treatment of precisely this model was given in 
[12]. In the following, we shall briefly summarize the previous results relevant for our 
considerations. 

The deterministic model, obtained by dropping the noise terms from (8), is given 
by the rate equations dtXi = at with a as in (12). Due to the normalization condition 
a + b + c = 1, its phase space can be viewed as the 2-simplex, i.e. an equilateral 
triangle. Its corners correspond to complete extinction (i.e. only one species of the 
three species is present), and its edges to states where one species is extinct and two 
are still present. The dynamics of the rate equations yields oscillations along closed, 
periodic orbits around a coexistence fixed point ata = 6 = c= |. Close to the fixed 
point, the orbits are almost circular, whereas further away, they approach the triangular 
shape of the simplex boundaries (see figure 3). These orbits are neutrally stable due to 



a 



( a{b 
b{c 
y c(a 



a) 

b) ) 



(12) 




(13) 
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the existence of a conserved quantity 



p — aba. 



(14) 



p assumes its maximum value pmax = 27 coexistence fixed point in the center of 

the phase space triangle, and its minimum value p = on the edges, corresponding to 
extinction of at least one species. 

In [12], it was also shown that with the inclusion of noise in the full stochastic 
model (4), the evolution of the population no longer takes place deterministically along 
one closed orbit, but can fluctuate randomly between different orbits (see figure 3). 
By means of a linearization around the coexistence fixed point, it was derived that 
eventually, the stochastic trajectory will hit one of the simplex boundaries, and from 
there move towards one of the absorbing corners of the triangle. This process means 
that two of the three species go extinct when stochasticity is included. It was also 
motivated that the mean extinction time scales as T oa N. 

Here, instead of linearizing the stochastic model (4), we will perform a stochastic 
averaging procedure over the deterministic orbits. This will remove the fast, oscillatory 
degrees of freedom (taking into account all non-linearities and the precise geometry of the 
phase space) and produce an effective one-dimensional stochastic differential equation 
for p. Through this, we will obtain an exact description of the extinction process and 
quantitatively correct results for mean extinction times. 

3.2. Stochastic Averaging 

Let us start with the formulation of the stochastic model using the Ito stochastic 
differential equations (8). 

Since the deterministic drift terms in (8) keep p = abc conserved, this quantity 
changes only due to the noise terms oc i.e. much more slowly than the oscillations 
along an orbit with constant p. Furthermore, p is a measure for closeness to extinction 
in the sense that the time when p becomes for the first time is exactly the time when 
the first of the three species goes extinct. Thus, looking at p allows us to separate the 
deterministic dynamics (i.e. the rapid oscillations along the closed orbits), which does 
not contribute to extinction, from the stochastic fiuctuations which lead to movement 
between different orbits and ultimately cause one of the species to die out. 

To determine the dynamics of p quantitatively, we can use (8) and apply the Ito 
chain rule [25], giving: 



The first term J^iOadiP is zero since p is conserved by the rate equations. (15) then 
implies that p changes on a slow time scale oc A^, and holds for any quantity conserved by 
the rate equations (9) of a stochastic Lotka- Volterra model. This shows in general that 
conserved quantities lead to neutrally stable coexistence in the classification of section 2. 
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The second term in (15) is a "stochastic drift" arising from the fluctuations and evaluates 
to: 

X 3 3 
— B^J^^^Jp = -—p. (16) 

The last term in (15) can be simplified by recalling that the sum of Gaussian variables 
is again a Gaussian random variable, with summed mean and variance. Applying this 
to the variables dW^ we can write J^k'^k dWk = o" dV with (^"^ = J2k'^k- ^ ^ ^'^^ ^ single 
Brownian motion with zero mean and unit variance. Inserting this into the last term in 
(15) yields: 

^ Y: {Cijd.p) dW, = -^a dV, (17) 

with = (Vp)C^C(Vp) = (Vp)5(Vp). 

Inserting this into (15) gives after some algebra: 
3 1 

dp^ -—pdt+ ^VOdV, (18) 

with 

D{a,b,c) = p^(-9 + - + l + -) . (19) 
V aba/ 

Note that (18) alone is not sufficient in order to describe the dynamics of p, since 
the coefficient D depends on all three variables a, b, c individually. The central idea 
of stochastic averaging [22] is now to replace the coefficient D{a,b,c) by the time- 
average -D(p) over the closed deterministic orbit with fixed p. This produces an effective 
stochastic differential equation purely in terms of p: 

dp=-^pdt+-^^/V{pjdV, (20) 

with 

Dip) := ^ j^^'^ D{a{t),b{t),c{t))dt. (21) 

Here T[p) is the period of a closed orbit of the rate equations corresponding to the value 
of p, and a(t), fe(t), c{t) parametrizes this orbit. D{a,b,c) is given by (19). 

Note that while superficially (20) looks quite similar to (18), their significance is 
quite different. While (18) describes the evolution of p as a function of the physical 
variables a{t), b{t), c{t), (20) is a closed equation for p as a free variable which we take 
to evolve on the interval 0; pmax = ^ • 

Intuitively, this procedure is justified by the time-scale separation described above: 
The deterministic drift along the orbit with constant p takes place on a time scale of 
0{1) fixed by the rate equations, and the movement between different orbits - i.e. the 
changes in p due to the noise terms in the stochastic differential equations (8) - occur 
on a time scale of 0{N). Hence, it can be assumed that the drift and noise terms for the 
evolution of p are in some sense "averaged" along the whole orbit before a significant 
change in p occurs. 
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0.002 




0.04 



Figure 4. Effective diffusion coefficient D{p) as given in (23). Observe that 
D{p) vanishes both at p = (extinction boundary of phase space triangle) and at 
P = Pmax = ^ (coexistence fixed point in the center of the phase space triangle). 



(20) can be reformulated equivalently as a Fokker-Planck equation for the 
probability distribution P{p,t): 

1 r„ 1 



d,[3pP{p,t)] + -d' D{p)P{p,t) 



(22) 



In the form (22), it is most apparent that the time scale of the extinction process is 
t (X N. 

In this simple model with its high amount of symmetry, the integral (21) can be 
performed analytically to give a closed expression for D{p) in terms of p. A sketch of 
the computation is given in appendix Appendix A, and the resulting formula is: 



D{p) = 3p' 



-3 + — + 
ai 



1 \ E{k) 



(23) 



/ K{k)_ 

Here, K{k) and E{k) are complete elliptic integrals of the first and second kind, 
respectively. The elliptic modulus is given by 



'^min)^l 



('^1 ^'min)'^max 

and amin, ctmax and ai arc the three real roots of the polynomial 
all — a)^ 



(24) 



4p. (25) 

As is well-known, these roots can be written down explicitely in terms of p. Graphically, 
D[p) is shown in figure 4. 

In total, with (20) and (22) we have provided a description of the extinction process 
in the rock-paper-scissors game as an effective one-dimensional stochastic process on the 
space of the deterministically conserved quantity p. This process has a linear drift term 
—3p and a complicated multiplicative noise given by (23). For comparison, we will now 
replace this by simple additive noise and verify the impact on the extinction process. 
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3.3. Constant Noise Approximation 

Replacing the complicated function -D(p) in (23) by a constant Dq, the extinction process 
(20) is reduced to a standard Ornstein-Uhlenheck process: 

dp=-^pdt+-^DodV. (26) 

In this approximation, the dependence of the mean extinction time on the starting value 
of p can be computed analytically in terms of the generalized hypergeometric function 

C is a constant fixed by the appropriate boundary conditions 

The full extinction probability distribution depending on time, starting from a 
fixed initial condition pq, can also be written down explicitely if the boundary at pmax 
is neglected. It then reads: 



PeAt, Po) = erfc l^poy ^^et _ i) j ■ (28) 

Asymptotically, this probability distribution possesses an exponential tail, Psurv(^) 
e~^*, independent of the precise value of Dq. The exponent —3t coincides very well with 
previous numerical results obtained in [30]. 



3.4- Comparison to simulations 

To verify the accuracy of the stochastic averaging procedure and the precise noise 
structure in (20), we simulated the mean extinction times and the extinction probability 
distribution of the original reaction system (10) using the Gillespie algorithm. The 
results are shown as blue dots in figures 5(a) and 5(b). 

This can then be compared to the predictions of the effective Fokker-Planck 
equation (22) with the full form of D(p) in (21) obtained by stochastic averaging. 
Although this effective Fokker-Planck equation cannot be solved analytically, the mean 
extinction times and survival probabilities can easily be determined numerically from 
the corresponding backward Fokker-Planck equation [25]. The results are shown as blue 
lines in figures 5(a) and 5(b). One can observe that the agreement to the simulation of 
the original reaction system (10) is excellent. 

Furthermore, figures 5(a) and 5(b) show the purely analytical results of the constant 
noise approximation in the previous section (namely, equations (27) and (28)) for 

^ From the singular nature of the boundary at Pmax = ^ in (22) and (23), one can derive that the 
mean extinction time must satisfy the boundary condition T'(pi„ax) = —9. In the constant noise 
approximation, this fixes the constant C to be: 
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(a) (b) 

Figure 5. Comparison of tlieory and simulation results for the rock-paper-scissors 
game. 

(a) : Mean extinction times depending on the initial condition. Solid curve: 
prediction obtained numerically from the stochastic averaging result (20). Dashed 
curve: Constant noise approximation (27). Crosses: results of direct simulation of the 
reaction system (10) using the Gillespie algorithm for N — 1500, averaged over 10'' 
realizations. 

(b) : Extinction probability distribution, starting from the coexistence fixed point 
p = dX t = Q. Solid curve: prediction obtained numerically from the stochastic 
averaging result (20). Red dashed curve: Constant noise approximation (28). Black 
dashed curve: Phenomenological approximation previously proposed in [12]. Crosses: 
results of direct simulation of the reaction system (10) using the Gillespie algorithm 
for N — 3000, averaged over lO'' realizations. 

Dq = 0.001. It can be seen that it reproduces the extinction probabihty distribution 
correctly and works well for the mean extinction times far away from the boundary at 
p = 0. Close to the boundary, the quantitative agreement is lost (which is not surprising, 
considering the shape of D{p) close to p = 0). This shows that close to the boundary, 
the precise form of the multiplicative noise plays a significant role, and provides further 
evidence for the correctness of (21). 

Having given an extensive treatment of the cyclic three-species model, we will now 
increase the number of species by one and consider the four-species model with cyclic 
dominance. 

4. Cyclic Four-Species Model 

In this section, we shall apply the formalism developed above to another example. We 
will consider the cyclic four-species Lotka-Volterra model, which is a natural object to 
study after the three-species rock-paper-scissors game discussed in the previous section. 
The reaction equations are given by: 

A + B ^ A + A, 
B + C ^ B + B, 
C + D ^ C + C, 
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D 



Figure 6. Phase space of the cychc four-species model. Green: Line of coexistence 
fixed points, (33). Red: Lines of extinction fixed points (34), (35), where only two 
non-interacting species remain. Grey: Closed deterministic orbits for various values of 
Ti, T2. Blue: Sample stochastic trajectory for N = 300. 



D + A ^ D + D. (29) 

For the case of equal reaction rates, the drift and noise terms for the Fokker-Planck 
equation (4) are obtained from (5) and (6) as: 

/ a{b -d)\ 



a 



B 



b{c — a) 
c{d - b) 
\ d{a - c) I 

( a{b + d) —ab 
—ab b{c + a) 



(30) 





\ —ad 



-be 



-be c{d -\- b) 







—ad \ 


—cd 

—cd d{a + c) ) 



(31) 



4.1. Rate Equations 

The rate equations for this model can be written down directly from the drift term (30) 
and read: 

/ dta \ ( a{b — d) \ 
b{c — a) 



dtXi = ai ^ 



dtb 
dtc 
V dtd I 



c{d - b) 
\ d{a — c) j 



(32) 
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As expected, the equations (32) keep the normahzation condition a + h + c-\- d — 1 
invariant. The phase space {a, 6, c, (i|a, 6, c, c? < Q;a + b + c + d— l}is now a 
three-dimensional simplex (i.e. a regular tetrahedron). Again, vertices correspond to 
extinction of all but one species, edges correspond to extinction of two out of the four 
species, and faces to extinction of one species. 

The fixed points of the rate equations (32) form three fines in the phase space 
simplex: 

• One line of coexistence fixed points given hy h = d and a = c, i.e. parametrized by 

r 



"•2 



(33) 



(a, 6, c, d) — ^i, - — t, - — with t e 

• The edge AC, i.e. all states with coexistence of A and C only, parametrized by 

(a, 6, c, (i) = (t, 0, 1 - t, 0) with t G [0; 1] . (34) 

• The edge BD, i.e. all states with coexistence of B and D only, parametrized by 

(a, 6, c, rf) = (0, t, 0, 1 - t) with te[0;l]. (35) 

A graphical representation of the structure of the fixed points is shown in figure 6. 

It is straightforward to check that the trajectories solving (32) now exhibit two 
conserved quantities: 

Ti = ac, 

T2 = hd. (36) 

The curves given by ri = const., T2 = const, are closed, neutrally stable orbits 
around the line of coexistence fixed points. Close to this line, they are almost circular, 
while further away they approach the shape of the simplex boundaries. A few exemplary 
orbits are shown in figure 6. 

In total, the deterministic dynamics for the four-species cyclic Lotka- Volterra model 
is quite similar to the behaviour in the three-species case. The relative abundances of the 
species oscillate indefinitely along a fixed, closed orbit in phase space, and no extinction 
takes place. 



4.2. Stochastic Extinction Process 

We would now like to investigate the behaviour of the four-species model when stochastic 
fluctuations, modelled by the noise terms in (8), are included. According to the 
experience from the three-species rock-paper-scissors model, we again expect to see 
extinction on a time scale oc N since the deterministic orbits are neutrally stable. 

However, as already apparent from the description of the rate equation dynamics in 
the previous section, the set of fixed points is now much larger than in the rock-paper- 
scissors game. Plugging the parametrizations of the fixed lines (33), (34), (35) into (31), 
we see that the noise matrix for the Fokker-Planck equation vanishes on the edges AC 
and BD, but not on the line of coexistence fixed points. 
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Hence, the absorbing states for the stochastic process are precisely the edges AC 
and BD of the phase space simplex, parametrized by (34) and (35), corresponding to 
states with a mixture of non-interacting species (A and C or B and D). This is consistent 
with the picture obtained directly from the reaction equations (29), where it is clear that 
these are exactly the states where no more reactions can occur. 

In order to analyze the stochastic model quantitatively, we now pursue the same 
approach used for the rock-paper-scissors model and investigate the behaviour of the 
deterministically conserved quantities when stochasticity is included. Applying the 
Ito formula and using the Langevin equations (8), we obtain the following stochastic 
differential equation for the conserved quantities Tq, a = 1,2: 

dTo.^^Y.iC^Ara)dW,. (37) 

Note that in contrast to the corresponding calculation in (18), there arc no stochastic 
drift terms due to the specific form of B in (31). Since {dri,dT2) is a two-dimensional 
Gaussian variable with zero mean (it is a linear combination of the Brownian motions 
dWj), it is uniquely determined by its correlation matrix B^. This is computed from 
(37) as follows: 

(dndn) (dndr^) \ ^ ( {Vn)B{Vn) {Vn)B{Vr2) \ 
{dr^dn) {dT2dT2) J [ (Vr2)5(Vri) (Vr2)5(Vr2) ) 

^ ( Ti{a + c){b + d) -4nT2 \ 

{ -4rir2 r2ia + c){b + d) J' ^ ^ 

Thus, we may equivalently write (37) as 

dra^4^j:c:^dV0, (39) 



B^ 



with chosen such that (C'^)^C'^ = B'^ given in (38). 

Note how this calculation provides a natural generalization of the analysis 
performed for the rock-paper-scissors game (section 3.2). As is apparent from (39) 
(or its Fokkcr-Planck equivalent), we again obtain a dynamics for the variables r which 
occurs on a time scale oc N. 

We furthermore observe that just as for the noise in the 3-species case, (38) still 
depends on the individual species concentrations a, b, c, d which vary along a trajectory 
with fixed ti, T2. To obtain a closed system of equations describing the evolution of ti, 
T2 we once again use stochastic averaging. 

As before, the high degree of symmetry and the simplicity of the model allow one 
to give an analytic expression for the average of the correlation matrix (38): 

B^K,.,)=f ^'/(-'•'■^) -f'-^ ), (40) 

\ -4TiT2 T2/(Ti,T2) J 



with / defined as: 



E{k{n,T2)) 



f(n,T2) := (a + c)(b + d) = MTi,r2) " ^ (41) 

K{k{Ti,T2)) 
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Again, K{k) and E{k) are complete elliptic integrals of the first and second kind, 
respectively. The elliptic modulus and the helper function h are given by: 



\ (7? — (To 

^Vi'^2)=4^-^, (42) 



1 



Ts) = - Ui + JfJi^ - (72 , (43) 



with 



(71 = 1 - 4ti - 4r2 (72 = 64Tir2. (44) 

For details of the calculation leading to these expressions, see appendix Appendix B. 

In total, we have now obtained a description of the extinction process in the four- 
species model of cyclic dominance by means of a two-dimensional stochastic process in 
the variables Ti and T2: 

1 2 

with T> chosen such that V^V — B'^ given in (40). As previously, ri and T2 are now 
treated as free variables. The region of ri-T2 space on which the process (45) occurs is 
bounded by the absorbing boundaries ri = and T2 — 0, as well as a reflecting boundary 
at 

v^+v^=^- (46) 

Equivalently to (45), we can write the effective stochastic process as a Fokker-Planck 
equation: 

dtP{T,, T2, ^) = ^ E ^cdf^ [(B^Unri. T2, t)] . (47) 

a, 13=1 

This provides a complete description of the extinction process in the four-species 
cyclic model as a two-dimensional diffusion with varying diffusion coefficients. 



4..3. Comparison to simulations 

As for the rock-paper-scissors game, we can now verify the accuracy with which various 
quantities of interest for the extinction process can be predicted by the effective 
stochastic process (45). Since this is now a two-dimensional stochastic process in a 
region with a complicated shape and mixed boundary conditions, it is much harder to 
treat than the one-dimensional effective process in the three-species case. 

Determining the mean extinction times and extinction probabilities from the 
effective Fokker-Planck equation (45) as was done in the three-species case is not feasible 
here, since it would require solving elliptic second-order PDE's over a domain bounded 
by (46). Instead, wc obtained mean extinction times and extinction probabilities from 
stochastic simulations of the cfTcctivc Langevin equations (45) using the XmdS package 
[31]. 
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Figure 7. Comparison of theory and simulation results for the four-species model of 
cyclic dominance. 

(a) : Extinction probability distribution, starting from a = b = c = d= j at 
t = 0. Solid curve: Simulation of effective stochastic process (45), averaged over 10* 
realizations. Crosses: Gillespie simulation of the reaction system (29) with a system 
size N = 8000, averaged over 10^ realizations. 

(b) : Mean extinction times depending on initial conditions. Lines: Simulations 
of effective stochastic process (45), averaged over 10* realizations. Crosses: Gillespie 
simulations of the reaction system (29) with a system size A'' = 2000, averaged over 
10"^ realizations. 



The results are shown in figures 7(a) and 7(b). It can be observed that the 
predictions of the effective Langevin equations (45) compare very well to the results 
of direct simulation of the original reaction system (29). We attribute the slight 
discrepancies in mean extinction times close to the boundary (46) to the general difficulty 
of simulating a stochastic process near a curved reflecting boundary. The stochastic 
averaging method becomes exact at this boundary, since the deterministic orbits are 
the individual coexistence fixed points. 

5. Conclusion 

In the previous sections, we have analyzed the extinction process in two stochastic 
Lotka- Volterra models which are neutrally stable in the deterministic treatment. We 
have seen that when fiuctuations arc included, the dctcrministically conserved quantities 
change slowly (on a time scale oc A^) and drive extinction. After some finite time, only 
non-interacting species remain. 

The separation of time scales between the rapid oscillations described by the 
deterministic rate equations and the slow movement between different orbits due to noise 
allowed us to apply the method of stochastic averaging. By doing this, we removed the 
fast, oscillatory degrees of freedom and gave a quantitative description of the extinction 
process using effective stochastic differential equations on the space of deterministically 
conserved quantities. 

We have obtained various quantities of interest for the extinction process from these 
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effective equations, and observed tfiat tfiey agree very well with direct simulations. 

The stochastic averaging procedure required computation of certain integrals over 
the closed deterministic orbits, which we were able to perform analytically in the toy 
models we considered. In more complicated models with less symmetry (e.g. different 
reaction rates), this may not be possible anymore. However, even for complicated and 
asymmetric closed deterministic orbits, the averaging required to determine the effective 
drift and diffusion coefficients can easily be performed numerically. Thus, we think that 
our approach should be just as helpful for elucidating the impact of noise in more general 
models. 

As we saw, a considerable advantage of the stochastic averaging method (especially 
in comparison to the treatment in [12]) is that neither the drift nor the noise terms 
need to be linearized. The full, nonlinear dynamics of the model and the multiplicative 
noise structure, as well as the complex geometry of the phase space can be taken into 
account. Furthermore, it is not necessary to write the dynamics explicitely in terms of 
a radial and a phase variable. This is useful since e.g. in the rock-paper-scissors model 
in section 3, there is no obvious choice for a canonical phase variable. 

The idea of describing the extinction process by the evolution of a deterministically 
conserved variable was also utilized in [32]. They apphed it in a semiclassical 
approximation and obtained the asymptotics of the extinction probability distribution 
in the standard two-species Lotka- Volterra model. However, our approach is quite 
different technically and allows a straightforward generalization to more complex models 
containing more than two species (as in the example in section 4). 

We also expect that it should be possible to extend this treatment away from the 
borderhne case of neutral stabihty to weakly stable or unstable models. Heuristically, 
this would give rise to a deterministic drift term in equation (20) (or its analogues) 
which is independent of N but controlled by some other small expansion parameter. 
Investigating how such models can be constructed in a natural way and the details of 
this generalization requires further research. 

Another approach to investigating the effects of stochasticity on similar models was 
pursued in [33, 34, 35], where the effects of the stochasticity on the phase variable and 
the spectral distribution of its oscillation were investigated. This is complementary to 
our treatment, where we focus on the radial variables instead. This allows us to capture 
the dynamics of the extinction process. 

In a more general context, the present discussion gives an illustration of how 
stochasticity may change the behaviour in a nonlinear dynamical system qualitatively 
by adding a stochastic drift to a deterministically conserved quantity. It also provides 
some ideas for treating the effects of complicated, multiplicative noise on such systems 
analytically. This may be of considerable interest for non-equilibrium statistical physics 
in general. 
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Appendix A. Computation of noise term for three-species cyclic model 

In this appendix, we shall sketch the computation of the stochastic averaging integral 
in (21). 

We will parametrize an orbit with fixed p by one of the species' concentrations, e.g. 
a. As is easily verified using the method of Lagrange multipliers, the extremal values 
Omin and Omax which a assumes on such an orbit are real roots of the polynomial 

a(l-a)2=4p. (A.l) 

The third root of this polynomial is then also real, and will be denoted by ai. As is well 
known, explicit expressions for all three roots exist. 
We hence write down the factorization 

a(l - a)^ - 4p = (a - amin)(a - ainax)(a - ai)- (A. 2) 

Now, let us give an explicit parametrization of each orbit. We will choose p and a 
as the independent variables, with < p < ^ and amin < a < Omax- Then b and c are 
given by 



a(l - a) ± ^a2(l - a)2 - 4ap 
= 2a ^^-'^ 

a(l - a) T ^a2(l - a)^ - Aap 
= L . (A.4) 

In each case, the + respectively - signs correspond to the two branches of the orbit for 
a fixed value of a. 

Plugging (A. 3) and (A.4) into the rate equation dta = a{b — c) we get: 

dta = ±^a2(l-a)2-4ap. (A.5) 
Now, let us calculate the period of an orbit, T{p). By a simple substitution we have 

T{p)= \dt = 2 ^. (A.6) 

The factor 2 arises since each orbit has two symmetric branches, when parametrized by 
e.g. a. Inserting (A.5) and (A. 2), we get 

T{p) = 2 / , (A.7) 

-'amin ^ a{a - amin) (a - Omax) {a - tti) 

This is a standard integral that can be expressed in terms of K{k), the complete 
elliptic integral of the first kind (see e.g. [36]): 
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The elliptic modulus k is given by (24). 

The last remaining piece we need is the average of ^ over a deterministic orbit. 
Again applying a substitution and using (A. 5) and (A. 2), we have: 

/■^(^) dt ^ /""max da ^ /■'^max da , . 

This, too, is a standard integral that can be expressed in terms of complete elliptic 
integrals (see e.g. [36]), giving: 

rnp) dt A [{a i-am^Mk^,k) + a^^K{k)] 

^max (^1 C'min)^max 



M w dt _ 
Jo a 



Here, E{k) is the complete elliptic integral of the second kind, Il{n, k) is the complete 
elliptic integral of the third kind, and the elliptic modulus k is again given by (24). The 
second line in (A.IO) was obtained by applying the relation T{{k'^,k) = 
Combining (A.IO) and (A. 8) we get the result used in (23). 

Appendix B. Computation of noise term for four-species cyclic model 

The computation of the average of the correlation matrix for the four-species model, 
(40), works along the same lines as the three-species case in appendix Appendix A. 

We parametrize the deterministic orbit with fixed ri, T2 in terms of a. The extremal 
values of a are given as: 

amax,min = ^ (l " 2 ± ^(1 " 2V^)' " 4Ti) . (B.l) 

The other variables are expressed in terms of a, ri and T2 as: 




1 - 

a 

c = -. (B.2) 
d 

As for the orbits seen in the three-species case, the plus and minus signs correspond 
to the two (symmetrical) branches of an orbit for each a. Prom (B.2) we get: 



dta = a(b - d) = ±^J{a{l - a.) - - Aa^T2. (B.3) 
With all of these results, we can now calculate the period of an orbit with fixed ri, 

T2: 

/■(Imax da 

T{n,T2) = 2 / , (B.4) 
^a„i, , / — a) — Ti)^ — 4a2r2 
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Again, the factor 2 arises from the two symmetric branches of each orbit. The 
denominator of the integrand can now be factored as: 

(a(l - a) - Ti)2 - 4a2r2 = ^(omax - a) (a - a^i^){ai -a){a- 02), (B.5) 

where 

1 



Note that ai > Omax and 02 < amin- Furthermore, Vieta's theorem gives following 
relations between Omin, Omax, 0.1, 02: 

aia2 = OtminO^max = 'Tl-i 

ai + 02 = 1 + 20^, 

amin + Omax = 1 - 2V^- (B.7) 

These will be very useful for simplifying some expressions later on. 

Plugging (B.5) into the integral (B.4), we get a standard elliptic integral 

r(ri,r2) = / K{k{n,T2)) . (B.8) 

Here, h is a helper function defined by 

h{Ti,T2) := (ai - ai„in)(amax - 02). (B.9) 
K{k) is the complete elliptic integral of the first kind, with the elliptic modulus given 

by 

e = '"7^^"^'^^^ . (B.IO) 

[ai — aniin)(«max ~ O2) 

By using (B.7), (B.l) and (B.6) we can express the helper function and the eUiptic 
modulus explicitly in terms of Ti and T2, giving the expressions (42) and (43). Observe 
that in (42) and (43), the symmetry in ti and T2 (which is required by the cyclic 
symmetry of the model, but was lost when we chose to parametrize the orbit explicitely 
using the variable a) is again manifest. 

To calculate the average of the noise matrix i?^ we also need the following integral: 



^a + c){b + d)^- la + -\(l-a--) dt. (B.ll) 

This is reduced to the following four basic eUiptic integrals: 

r^max a da 



ran 
ran 

J Ami: 

ran 

J ami 



h 
h 



^(ai - 


a)(amax - a)(a - 
da 


O^min) ifl 


a2) 


a\J{ai - 


- a)(amax - a)(a - 
0^ da 


~ ^min) ifl 


- 02)' 


sjiai - 


a)(amax - a)(a - 
da 


^^min) ifl 


02)' 


a2^(ai 


- a)(amax - a){a 


Q'min) ifl 


- 02) 
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These can be computed using the formulae in [36] giving: 

2 r 

/i = ^ aiK{k) - (ai - amax)^(Q^^ k) 



h 
h 

/4 = 



-K{k)-{--^]Ii{ai,k) 



y/h 
2 

2 r 

-j= alK{k) - 2ai(ai - a^^)U{a'^, k) + (ai - Omax) ^2(0;^ 
2 



1 2/1 1^^ 

Vh 

As usual, K{k), E{k) and n(n, fc) are the complete elliptic integrals of the first, second 
and third kinds, respectively, a, cti and V2 are given in our notation by: 

o ^2'Tnax ^min 



a — , 



^1 1 O-max 



«1 - 7 Z V- 

V2 is defined as 

^^^"^^ 2(a2-i)(P_«2) ■ (B-12) 

Here, we dropped the elliptic modulus k (which is always the same) from the arguments 
of K, E and 11. Combining the expressions for Ii,...,!^, we obtain after some long 
and tedious algebra the surprisingly simple result (41) for (B.ll). This result, too, is 
symmetric in ri and T2 (as expected, since the quantity {a + c){b + d) is invariant under 

cyclic permutations). 
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